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Abstract 

We propose a novel model reduction approach for the approximation of non linear hyperbolic equations 
in the scalar and the system cases. The approach relies on an offline computation of a dictionary of 
solutions together with an online Zd-norm minimization of the residual. It is shown why this is a natural 
framework for hyperbolic problems and tested on nonlinear problems such as Burgers’ equation and the 
one-dimensional Euler equations involving shocks and discontinuities. Efficient algorithms are presented 
for the computation of the ZJ-norm minimizer, both in the cases of linear and nonlinear residuals. 
Results indicate that the method has the potential of being accurate when involving only very few 
modes, generating physically acceptable, oscillation-free, solutions. 


1 Introduction 


Many engineering applications require the ability to simulate the behavior of a physical system in real-time. 
This requirement holds in particular when a full parametric exploration of the behavior of the system is 
sought. In aerodynamics, such an exploration can be done to compute the flow around an aircraft for varying 
boundary conditions or to design its shape to maximize lift and minimize drag. Uncertainty quantification 
also requires a large number of simulations with varying parameters in order to propagate chaos by means 
of a Monte-Carlo method or calibrating input parameters by a Markov chain technique. A third important 
application is flow control. 

When such a large number of simulations is required, the cost of one simulation is critical to the ap¬ 
plication at hand. This cost can be lowered by using sophisticated computer science techniques such as 
parallelization but such techniques are usually not enough to allow full parametric exploration, especially 
when computational resources are limited. 

Alternatively, model reduction techniques can alleviate the cost of such repeated simulations with limited 
computational resources Bill- Model reduction is directly based on the underlying high-dimensional 
model (HDM) that results from a standard finite element, finite volume of finite differences formulation. In 
the present paper, Partial Differential Equations (PDE) of the following type are considered: 


dU 

~dt 


+ L(U) = 0 


B(U) = g 


U(x, t = 0) = Uq(x ) 
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ieO, t e [o, t\ 
xedn, t e [o, t] 
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(1) 


L is a differential operator (for example the Laplacian or the divergence of a flux), and B a boundary 
operator, fn this paper, we are particularly interested in the case where the solution U(x,t) £ is a scalar 
or a vector and L is the divergence of a flux F. Two examples will be considered by increasing order of 
complexity: 

• Burgers’ equation for which U = u is scalar: 


— Its unsteady version, 

du d (1 2 \ 

» + & U"') = “ (x ' 0)= "” w 

with periodic boundary conditions 

— It steady version with weak Dirichlet boundary conditions 

• The one-dimensional compressible Euler equations for which U = (p, pu, E ), F(U) = ( pu , pu 2 +p, u(E+ 
p)) and the perfect gas equation of state holds: 

P = (7 - !) (e - \pu 2 ^j ■ 


p denotes the density, u the velocity, p the pressure and E the energy. 

After discretization in space, the solution is denoted as u(t) £ WL Np . The PDE is here parameterized by a 
parameter vector /x £ R m that allows changes in the operator L , the boundary operator B or the initial 
conditions. For simplicity and without loss of generality, this parametric dependency will be omitted in the 
next paragraphs. 

Instead of allowing any value of the solution degrees of freedom u, model reduction however restricts the 
solution to be contained in a subspace of the underlying high-dimensional space. This subspace is determined 
by an optimized reduced basis that is determined in a training phase. Thus, a large number of degrees of 
freedom (say millions) are represented by only a few number of coefficients in the representation of the full 
solution in terms of the reduced basis vectors, leading to important computational savings. Two important 
questions arise at this point: (1) How can an optimal reduced basis be constructed? and (2) How can the 
evolution of the reduced coefficients be computed in a stable fashion? 

A popular method for choosing an “optimal” basis is Proper Orthogonal Decomposition (POD), first 
introduced as a tool for the analysis of flows by Lumley [5] and then extended and popularized by Sirovich [5] . 
The idea behind POD is to collect a few snapshots of the solution and then compute the best approximation of 
these snapshots in terms of a small number of reduced basis vectors. Mathematically speaking, if Uj(f/) £ R p 
denotes the value of the discrete solution u at grid point x,, i = 1, ■ • • , N and at time tj, l = 1, • ■ • , Nt, 
POD constructs M orthogonal functions <p e £ [L 2 (R d )] p such that the following functional is minimized: 


N t Np 

J&ir -.&*) = ££ 


1=1 i=1 


M 


Ui(ti) - ) (u(ti),(t>t)<p ei 


1=1 


( 2 ) 


where <$>&£ R p denotes the value of <f> at x,;. || • || denotes here the Euclidean norm in R p , and ( • , • ) is the 
L 2 norm. A minimum of the functional J can be analytically computed by Singular Value Decomposition 
[T, and the reduced basis vectors are the left singular vectors of the snapshots matrix 


( Ui(fi) ... Ui(tjVt) 

Ujv(fi) ... u N (t Nt ) 


Defining by {A^}^ the positive eigenvalues of S T S sorted decreasingly, the error associated with the mini¬ 
mum of the functional is 

N t 

Wi,--- ,0 m)= E (3) 

e=M +1 
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In the continuous case, the functions </>^(x) £ R p , are the solution of Fredholm alternative 

f I?(x, x')</> ^ (x , )dx , = A^0^(x), for all x £ fl, (4) 

An 

where R(x,x') = u(x)u(x') T . 

In both the discrete and continuous cases, the basis dimension M is determined on the basis of the 
decay of the eigenvalues A^. Given a tolerance e < 1, M is selected as the smallest dimension such that the 
following relative truncation error is smaller than e. 



,<Pm) 

(5) 



In general, one expects the eigenvalues Xe to decrease very rapidly to 0. This allows, when this assumption 
is true, to consider only the most energetic modes in the decomposition. Unfortunately, it is not always the 
case that the eigenvalues Xe are rapidly converging to zero. This is demonstrated by the following simple 
counter example for which a simple scalar advection problem defined on f 1 = [0,1 [ is considered: 


du du 

dt dx 

(6a) 

with the boundary condition 

u{ 0, t) = 1 

(6b) 

and the initial condition 

u(x, 0) = 0. 

(6c) 


The solution is given by a traveling discontinuity 


f 1 if x < niin(t. 1) 

^ | 0 otherwise. 

Considering grids Xi = i/N, i = 0,..., N for varying number of grid points N and snapshots collected at 
times as tk = kAt, with At = 1/N, a series of POD bases is constructed numerically. For each grid size N, 
the eigenvalues A t(N) are reported in Figure [l] One can observe that the ratio Xe(N)/Xi(N)) behaves like 
1/fc and max(A^) behaves like 0.63 A. This illustrates that it is not possible to select only a few dominant 
modes, due to the slow decay of the POD eigenvalues. This example also illustrates why most of the work 
on model reduction has been focused on regular problems, and for fluids, on incompressible flows, see e.g. 
among many others (HI ©i, HU- For compressible (but regular) flows, one of the early work is [IT], then one 
may mention m for compressible turbulent flows, m for compressible inviscid flows and PI H51H51 13] for 
the case of linearized compressible inviscid flows. 

Concerning compressible fluids, there is another difficulty. In problem Q, one needs a norm. In the case 
of incompressible flows, a natural norm is related to the kinetic energy. For compressible materials, however, 
one needs to take into account the density, velocity and the energy, i.e. the thermodynamics. A simple 
T-norm cannot be used because one cannot combine in a quadratic manner these variables, for dimensional 
reasons. Only a non-dimensionalization of the variables m can alleviate the dimensionality issue mmm- 

The natural equivalent of the Z, 2 -norm is however related to the entropy, which is not quadratic: if a 
minimization problem can be set up, its solution is non trivial. These arguments were raised in m, and an 
energy-based norm was developed in mm for linearized compressible flows. 

To circumvent those issues, an approach based on a dictionary of solutions mm is developed in this 
work as an alternative to using a truncated reduced basis based on POD.The elements of this dictionary are 
solutions u (ti\/ij) computed for varying values of time ti and parameter n J £ Selecting appropriate 
parameter samples /ij £ T> C R m is a crucial step that can affect the accuracy of the reduced-order model in 
the parameter domain. Greedy sampling procedures have been developed when error estimates are known EB 
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Figure 1: Ratio of POD eigenvalues \og(\k{N) /Ai(iV)) for N = 400, 600, 800,1000,1500 grid, points. 

mmmmmm- In this work, the issue of optimal sampling is not considered the main focus of the paper 
lies in establishing an effective model reduction approach based on dictionaries for hyperbolic problems. The 
development of a strategy to sample optimal values of /x in this context will be the topic of further research. 

In addition to choosing an appropriate dictionary V, selecting an approach for computing a reduced 
solution based on that dictionary is also crucial. For self-adjoint systems, Galerkin projection is a natural 
approach but there is no motivation for using Galerkin projection for nonlinear compressible flows. Instead, 
strategies based on the minimization of the residual arising from the reduced approximation have been 
successfully developed for compressible flows in dJ n ED Eg. These approaches rely on a minimization of 
the residual in the L 2 sense. In the present work, this minimization problem is extended to the more general 
minimization using a L 9 -norm, with emphasis on q = 1 and q = 2. For nonlinear systems, an additional step, 
hyper-reduction, is required to ensure an efficient solution of the reduced system E51 dm3]. Hyper-reduction 
is not considered in this work but will be the subject of follow-up work. 

This paper is organized as follows. Motivations for using the Zd-norm in the case of hyperbolic systems 
are given in Section [2] where we show that q = 1 is very closely linked to the concept of weak solutions 
of hyperbolic problems. The proposed model reduction approach is then developed in Section [3] in both 
the steady and unsteady cases. Finally, the proposed procedure is applied to the model reduction of several 
steady and unsteady systems in Section [4] and conclusions are given in Section [5j 

2 Motivation for the L 1 -norm. 

In solving minimization problems, it is quite usual to minimize a residual with respect to the Z 9 -norm for 
a suitable q. The choice q = 2 is very common because it amounts to minimize in some least-squares sense 
and many efficient algorithms are available. In the case of hyperbolic problems, as we are concerned with 
here, this is still a convenient choice (after proper non-dimensionalization as mentioned above), but it might 
not be the most natural one, as demonstrated in the work of Guermond et al. on Hamilton Jacobi equations 
and transport problems (26[ l27| . In particular these works show, at least experimentally, that the numerical 
solution has an excellent non-oscillatory behavior by minimizing the ZJ-norm of the PDE residual. In fact, 
this observation is our original motivation for choosing the Z 1 -norm, since we are interested in preserving 
the non-oscillatory nature of solutions. In this section, we further justify the choice of the L 1 -norm applied 
to the residual, and show that it is closely related to the weak formulation of the problem. The following 
discussion is formal. 


4 







Let us consider the problem 


( 7 ) 


dU 

— + div F(U) = 0 

defined on C R d , t > 0. The steady problem can be done in the same exact manner. We assume that the 
solution U belongs to R p , so that F = (Fi ,..., F p ) T . The weak form of this is: for any <p £ [C 1 (f2)] P and 
with compact support, we have: 


J ip(x, t ) (Y + div F(U)j dx = 0. 


Integrating by parts yields 


[ Udx+ [ V<p ■ F(U)dx = 0. 

Jn ot J n 

If we restrict ourself to the set of test functions {ip £ [C' 1 (fI)] P , |M|oo < 1}, we have that U is a solution if: 


sup 

{^[CTnhMIvlU^i} 


f Udx+ f Vip ■ F(U)dx \ = 0. 
Jn Jq J 


Let us now real the definition of the total variation of a function g £ L 1 


TV(g) = sup 

v ec 0 Hii^n£~(R d ),IMIoo<i 


/ Vc p{x ) • g(x)dx l, 
J R<* J 


and we see that if in addition g £ C 1 (K d ), TV(g) = f Rd ||Vg||dx = || 

Before going further, let us mention the following classical result that will be useful. Consider {xi}i e % a 
strictly increasing sequence in R, we define x i+ i/ 2 = x,+x '+ 1 . We assume that R = Uigz[®*_ 1 / 2 ) a; *+i/ 2 [ an d 
consider g defined by: for any igZ, 


g(x) = gt if x £ [xi-x/ 2 ,x i+1 / 2 [, 

we see that _ 

TV id) = Y 1 9i+i ~ 9i\- 

iez 

Thanks to this definition, we see that if we define the space-time flux F = ( U,F ), U is a weak solution 
if and only if the total variation of F vanishes, TV (F) = 0. 

Now, instead of having the exact solution, we consider an approximation procedure that enables, from 
u n ss U{ . ,t n ), to compute u n+1 ss U( . ,tn+ 1 ), sa y £(u",u” +1 ). 

For instance, assume that we have a finite volume method and d = 1: for any grid point i £ {1, • • • , N}, 

[£(u n , u" +1 )]. = Az« +1 - <) + Af(f i+1/2 K) - fj- 1 / 2 (u")). 

A way to evaluate u" +1 is to minimize the total variation, i.e. 

TV(C) = Y |A^« +1 - O + At(f i+1/2 (u") - fi_ 1/2 (u"))|, 

iS Z 

A x{-Vi - u”) + At(f i+1/2 (u n ) - fi_ 1/2 (u")) . 

Clearly, if I is equal to the set of grid points, the solution is given by 

< +1 = u " - ^ (^fi+i/2(u n ) - f J -i/ 2 (u 11 )^ 


u " +1 = 


argrnin Y^ 

v piecewise constant 
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and nothing new is gained. When X is not equal to the set of degrees of freedom, then something new 
happens. We expect precisely to exploit this idea, or ideas related to this. 

In the remainder of this paper, this idea is exploited in the case of model reduction, for which X is not 
equal to the set of grid points and the TV semi-norm slightly modified in order to guarantee (1) that a 
unique solution to the minimization problem exists, and (2) that the minimization problem is as easy as 
possible to solve. 


3 Formulation 

3.1 High-dimensional model 

Without loss of generality, the case of the classical finite volume method is considered to define the High 
Dimensional Model (HDM). A computational domain il C K d is considered, and in most of this paper, 
Q C R, that is d = 1. Starting from a subdivision ... < Xj < Xj+\ < . .., we construct control volumes 
Kj = [xj — 1/2 j Xj -j-i /2 [, jeZ where 

x j —t - x j r 
Xj+l/2 = -X-• 


A finite volume semi-discrete formulation of 0 writes 

. du 


Ki\-£ + f i+!/2( u ) - f?-l/2( u ) - 0 


(8a) 


where fj- 1 - 1/2 is a consistent numerical flux. In each applications, we consider Roe’s formulation and a first 
order scheme. We assume either compactly supported initial conditions or initial conditions with periodicity 


u j(t = 0) 


\k 3 \ 


Uo(x)dx. 


(8b) 


/if. 


In (8a I, Uj stands for an approximation of the average of the solution in the cell Kj , 


Uj(i) 


1 

Wj\ 


U ( x , t)dx. 


IK, 


The time stepping is done in a standard way, for instant by Euler time stepping. 


3.2 Model Reduction by residual minimization over a dictionary 

3.2.1 Steady problems 

Two approaches are available to solve a steady state associated with problem (jT|) . The first one is to use a 
homotopy approach [E5| with pseudo-time stepping, resulting in the solution of an unsteady problem which 
limit solution is the desired steady state. The procedure described for unsteady systems in Section |3.2.2| can 
be, in principle applied to this case. The second approach is by a direct solution of the steady-state problem. 
The discretized steady-state problem writes 


r(u (/x),M) = 0 

where r(-, •) is usually a nonlinear function of its arguments, referred to as the residual. This set of nonlinear 
equations is typically solved by Newton-Raphson’s method. This second approach is followed in this work 
for steady problems. 

The parameter vector ji £ V C R m can, for instance, parameterize the boundary conditions associated 
with the steady-state problem. The parametric domain of interest V is assumed here to be a bounded set of 
R m . 


6 








The solution manifold M. = |u(^i) s.t pi £ V C K m } is assumed to be of small dimension. This manifold 
M belongs to L°°(R d ) n I?U(R d ), and thus can be locally described by some mapping 9 : V >->■ L°°{R d ) D 
BV (K fi ).To approximate this mapping, we consider a family of r parameters in V, {/r f }£ =1 , and compute 
the associated dictionary of solutions V = {u (fx e )} r (=1 of (| 8 |. 

The steady-state u(/x) is then approximated as a linear combination of the pre-computed dictionary 
elements T> as 

r 

u (m) ~ ( 9 ) 

e=i 

For a new value of the parameters fi £ V, the reduced coordinates {a^(/i )}^ =1 are then computed as the 
solution of the minimization problem 


a(/x) := (ai(Ai),...,< 


= argmin 

@={(h >■ 



,p 


( 10 ) 


In this paper we consider for J the following convex functionals, which are described in more details in 
Appendix [A] 

• the L 2 -norm J(z, x) = ||z || 2 and its regularized version J(r,/3) = ||r || 2 + 77 ||/ 3|| 2 with 77 > 0. 

When r is a linear function of (3 (r = A/3 + b), this choice of functional results in the solution of 
an ordinary least-squares problem, described in Appendix |A.1.1| When r is a nonlinear function of 
/3, Gauss-Newton or Levenberg-Marquardt procedures can be used to minimize J, as described in 
Appendix |A.1.2[ 

• the L 1 -norm J(r,/3) = ||r||i or its regularized variant, J(r,/3) = ||r||i +7/||/3||i with 77 > 0. 

Two approaches are considered to minimize J when r is a linear function of /3. 


1. Linear Programming (LP), involving the solution of an optimization problem with 2m+r variables 
and 3?n constraints. 

2. The Iteratively Reweighted Least-Squares approach (IRLS) [29] . 


Both approaches are described in great detail in Appendix | A.2. 1| When r is a nonlinear function of /3, 
a Gauss-Newton-like procedure can be used in combination with either the LP or IRLS approaches, as 
described in Appendix |A.2.2| Unicity of the solution can be guaranteed by setting the regularization 
term 77 > 0 . 


The Huber function J(r) = )C)=i [30] as described in Appendix A.3 In the present work, 


minimization of the Huber functional is carried out by the IRLS approach. The procedure is described 
in details in Appendix |A.3| 


The Huber functional can also be regularized by defining J(r,/3) = Y^=i ^M^i) + ??||/3 || 9 with q = 1 
or q = 2 and 77 > 0 . 


Remark 3.1. • Decreasing the dimensionality of the solution space from N to r is not enough to gain 

computational speedup when the system to be solved is nonlinear. An additional level of approximation, 
hyper-reduction, is necessary jumm. Hyper-reduction is not considered in the present work, 
but will be the focus of future work. 

• A careful selection of the sample parameter samples is necessary in order to generate a reduced- 

order model that is accurate in the entire parameter domain V. Greedy sampling techniques [F0 i \lfK 
0 2f[ \22l 12,91 S| 1 24V , associated with a posteriori error estimates, have been successfully used to construct 
reduced models that are robust and accurate in a parameter domain V. These techniques are not 
considered in this paper but will also be the focus of future work. 
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3.2.2 Unsteady problems 

For simplicity, in the remainder of this section, we assume that only the initial condition u°(/j) depends on 
a parameter vector fi £ V C Again, the family of solutions u(/x) of the Cauchy problem © is then 
conjectured to belong to a low dimensional manifold M. when the initial condition is parameterized in (8b I. 


To approximate this mapping, we consider a family of r parameters in V , {Mcj+i) and compute the 
associated solutions of (|8j) for respective initial conditions u°(p f ), t = 1,..., r. 

Once these solutions are computed, we propose to approximate, for any parameter fi £ V the solution 
{u"(M)}^o associated with an initial condition u°(/x) by approximating it as 


u”(+ = J2<x n {n)u n (n<>) 


l 


by the following procedure: 

1. Initialization: determine the reduced coefficients {a°(/x)}£ =1 as: 


£*°(m) := (+(/+■••,++)) = argmin V/WV+/3 

/3=(/3i,-,/3 r ) 

for a given choice of functional J(u,/3). 

2. Assume that a n (fi) = (a”(/x),..., a”(/x)) is known, determine a" +1 = (a™ +1 ,..., a” +1 ) such that: 




n+1 


(M) = ^ argmin . J ^EA u " +1 W - w n (fx) - ^ (Vi/ 2 (w n ) - f_ 1/2 (w n ) j 


where 


w "0) = 


We see that the second step can be written as: find a n+1 (fj,) that minimizes 


J{A n+1 a n+1 - b") := J(A 
where the matrix A” +1 can be written by blocks as 

K +1 (Mi) 


n+1 a n+1 - b n 


r i+l 


A” +1 = 




n+1 


n+1 

l N 


(P'rT 


(Mr)y 


( 11 ) 


and b" depends on known data. 

A few immediate remarks can be made. 

Remark 3.2. • In the case of a linear flux, Problem ([!]) is linear. If St. is the mapping between the 

initial condition Uq and the solution at time t, we have S t (U + V) = S t (U) + S t (V). This means the 
exact solution of the Cauchy problem with Uq = (/+ is SflUo) = a®St.(Uo{pii))■ In the 

case of a linear scheme, minimizing the functional J should result in a. n = cc 0 for any n > 0. 

• In the case of an explicit background scheme, the choice of the numerical flux, how high order is reached, 
and the choice of time stepping has no influence on the overall procedure: any sub-time step would be 
treated similarly. In this paper, we have chosen a first order method with Euler time stepping in the 
case of unsteady problem. 


In the case of an implicit scheme, a Newton-like procedure can be applied to minimize the functional as 
in m- At each time step, the procedure is then identical as in the steady case described in Section\3.2.1\ 





4 Numerical examples 

4.1 Simple regression 

As a first example, the choices of functionals proposed in Section [3.2.1| are applied to a very simple regression 
problem. This example illustrates the behavior of each approach. In a first case, 22 points {cE;}f£i are 
randomly selected in the interval [0,1] and the coordinates {y,} 2 '^ are drawn from a distribution 2 Xi + 
0.4 + 10 -1 U{— 1,1) where U{— 1,1) denotes the uniform distribution between —1 and 1. The regression 
approximation is therefore y ss a\x + a .2 and the target is a* = (a*, a^) = (2,0.4). 

The Id-norm minimization by LP and IRLS, L 2 -norm and Huber function minimizations are then used 
as functionals for that problem. The results are reported in Figure [2ja) and Table [T] One can observe that 
all four methods provide a good approximation of a*. Furthermore, the two L 1 minimization procedures as 
well as the Huber function minimization return identical results. 



Target 

L 2 -norm 

L 1 -norm (LP) 

I 1 -norm (IRLS) 

Huber function 

ai 

2 

1.9520 

1.9037 

1.9037 

1.9037 

a 2 

0.4 

0.4087 

0.4408 

0.4408 

0.4408 


Table 1: Regression: results without outliers 



Target 

L 2 -norm 

Z/ 1 -norm (LP) 

I 1 -norm (IRLS) 

Huber function 

ai 

2 

0.9256 

1.9037 

1.9037 

1.9037 

Q-2 

0.4 

0.9545 

0.4408 

0.4408 

0.4408 


Table 2: Regression: results with outliers 


In a following set of experiments two outliers points are defined and all four approaches applied to that 
new regression problem. The results are reported in Figure [2][b) and Table [2] One can observe that the 
L 2 -norm minimization procedure is much more sensitive to the outliers. As such, the regression coefficients 
returned by that procedure differ greatly from the previous case and a* is inaccurately estimated. On the 
other hand, the Id-norm and Huber minimization procedures are much less sensitive to the outlier points 
and accurate estimations of a* are provided. Again, the three estimations are identical. 

4.2 Model reduction of steady problems 

4.2.1 One-dimensional advection equation 

The following one-dimensional steady advection equation is considered: 

nil 

^(x) = f(x-,y),x&[ 0,1], (12) 

where 

,, . _ ~2k exp(—2k(x — y)) 

X '^ (l + exp(—2 k(x — y))) 2 

and k = 100. The solution exhibits a sharp gradient at location x = y similar to a shock. A Dirichlet 
boundary condition u{ 0) = 1 is imposed at x = 0. This PDE is discretized by finite differences using a 
uniform mesh, resulting in a HDM of dimension N — 10 3 . 

A dictionary of r = 6 solutions is built for y G {0.3,0.34,0.38,0.42,0.46,0.5}. The six solutions are 
depicted in Figure [3j Each solution has a high gradient at a different location x = y. Five different model 
reduction methods are then compared, namely Galerkin projection, I 1 -norm minimization by LP and IRLS, 
L 2 -norm and Huber function minimizations. 
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2.5 



(a) Case without outliers 



(b) Case with outliers 

Figure 2: Regression: Comparison of the methods 


10 
















Figure 3: One-dimensional advection equation: dictionary and target solution 



-Exact 

- Galerkin 

-L2 

-LI by Linear Programming 

LI by Reweighted Least-Square;; 
-Huber 


Figure 4: One-dimensional advection equation: solutions at one of the dictionary parameters n* = 0.34 
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Figure 5: One-dimensional advection equation: predicted solutions at target parameter fi* = 0.45 


In a first experiment, a target parameter fi* = 0.34 belonging to the dictionary is selected. As reported 
in Figure [4] all approaches correctly predict the solution. Next, a target parameter fi* = 0.45 not belonging 
to the dictionary is selected. The model reduction results are reported in Figure [5] One can observe that the 
L 1 -norm and Huber function minimizations lead to predictions that are identical and are the most physical 
as they exhibit a single shock. On the other hand, model reduction based on Galerkin projection predicts 
an unphysical solution with very large oscillations. L 2 -norm minimization leads to a smooth solution that is 
not physical either. The reduced coordinates corresponding to each of the six model reduction approaches 
are reported in Table [ 3 ] One can observe that the solutions returned by the approaches based on iP-norm 
and Huber function minimizations are very similar and are sparse, unlike the solutions returned by Galerkin 
projection and L 2 -norm minimization. Finally, the residuals returned by each minimization technique are 
depicted in Figure [6] One can observe that L 2 -norm minimization results in a higher maximum local residual 
as well as non-zero residuals that have a much larger support in the computational domain. 



Galerkin 

L 2 -norm 

T 1 -norm (LP) 

L'-norm (IRLS) 

Huber function 

Oi\ 

-0.896 

0.046 

-3.6xl0" 12 

-2.5x 10" 11 

-1.4xl0” lu 

012 

0.962 

0.045 

5.0x10"® 

4.3x10-® 

1.4x10"® 

03 

-1.028 

0.045 

-2.0xl0" 5 

-1.8xl0" 5 

-7.4X10" 6 

014 

1.115 

0.097 

0.033 

0.031 

0.019 

C*5 

0.417 

0.725 

0.967 

0.970 

0.981 

C *6 

-0.001 

0.040 

-4.6xl0" 4 

-5.0xl0" 4 

-2.9x 10" 4 


Table 3: One-dimensional advection equation: reduced solutions 


4.2.2 Two-dimensional advection-diffusion equation 

The two-dimensional advection-diffusion equation is then considered 

A (n) ■ Vu(:r, y) — kA u(x, y) = 0, (x, y) £ Cl = [0,0.018] x [0,0.018] (13) 
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Figure 6: One-dimensional advection equation: residuals at target parameter //* = 0.45 

with incoming Dirichlet boundary conditions and outgoing Neumann boundary conditions. The problem is 
parameterized by the angle of the advection flow with respect to the x axis: A(/i) = (||A|| 2 cos(/r), || A|| 2 sin(/r)). 
This problem is dominated by advection since ||A|| 2 = 0.5 and k = 2 x 10 - '. 

The problem is discretized by finite differences using a uniform mesh with 304 points in each direction, 
resulting in N = 88464 degrees of freedom. For this large scale problem, solving the L 1 -norm minimization 
problem by LP is not tractable and as such only the IRLS method is used in the Zfl-norm case. 

A dictionary V of two solutions is constructed for /r £ {f, f } and a target parameter /i* = ^ considered. 
The respective solutions are depicted in Figure [7] 

The following four model reduction methods are then applied: Galerkin projection, L 2 -norm minimiza¬ 
tion, L 1 -norm minimization by IRLS and Huber function minimization. The corresponding reduced solutions 
are reported in Table [4] and the solutions and errors in Figures [8] [ 9 ] For this problem, the L'-norm method 
by IRLS failed to converge and the returned solution is zero. However, the Huber function minimization 
approach was much more robust and returned a physical solution with sharp gradient. Hence, this example 
illustrates the advantage of using the Huber function versus pure L 1 -1101111 minimization. Galerkin projection 
and L 2 -norm minimization returned very similar but much less physical solutions with gradients that are 
much less sharp. 



Galerkin 

L 2 -norm 

L i -norm (IRLS) 

Huber function 

OL\ 

0.543 

0.467 

1.9xl0 -11 

0.021 

a 2 

0.688 

0.529 

5.2xl0 lu 

0.979 


Table 4: Two-dimensional advection-diffusion equation: reduced solutions 
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Figure 7: Two-dimensional advection-diffusion equation: dictionary and target solutions 
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(b) L 2 -norm minimization 

Figure 8: Two-dimensional advection-diffusion equation: predicted solutions and errors at target parameter 
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Figure 9: Two-dimensional advection-diffusion equation: 
(continued) 


predicted solutions and errors at target parameter 
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Figure 10: One-dimensional Burgers’ equation: dictionary and target solution 


4.2.3 Steady Burgers’ equation 

The following one-dimensional steady Burgers’ equation is considered: 


1 8{u 2 ) 

2 dx 


(x) = f{x-,n), ie!l = [0, 1], 


(14) 


where 


f{x-,n) = 


—2fcexp(— 2k(x — fi))( 1 + 3exp(—2 k(x — n)) 
(1 + exp(—2 k{x — n))Y 


with k = 100. The solution exhibits a strong gradient at location x = ji. A Dirichlet boundary condition 
u(0) = 1.5 is applied at x = 0. This PDE is discretized by finite differences using a uniform mesh, resulting 
in a HDM of dimension N = 10 3 . 

A dictionary of r = 6 solutions is built for /i £ {0.3,0.34,0.38,0.42,0.46,0.5} by solving each steady 
state problem by Newton-Raphton’s method. The six solutions are depicted in Figure [lOj Each solution 
has a “shock" at a different location. For this case, six different model reduction methods are compared, 
namely Galerkin projection, L 2 -norm minimization by Gauss-Newton and Levenberg-Marquardt, L-norm 
minimization by LP and IRLS and Huber function minimization. 

A target parameter /i* = 0.45 not belonging to the dictionary is selected. The model reduction results 
are reported in Figure 11 One can observe that the L 1 -norm and Huber function minimization results are 


identical and are the most physical as they exhibit a single “shock". On the other hand, Galerkin projection 
results in an unphysical solution with very large oscillations and inaccurate constant solutions before and 
after the “shock". The two L 2 -norm-based approaches result in a very smooth solution before the shock 
and an undershoot after the shock that are not physical either. Finally, the residuals returned by each 
minimization technique are shown in Figure 12 One can observe that L 2 -norm minimization approaches 


result in a higher maximum local residual as well as non-zero residuals that have a larger support in the 
computational domain. 


17 



















Figure 11: One-dimensional Burgers’ equation: predicted solutions at target parameter /j,* = 0.45 
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Figure 12: One-dimensional Burgers’ equation: residuals at target parameter fi* = 0.45 
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4.3 Model reduction of unsteady problems 

4.3.1 Unsteady Burgers’ equation 

We consider here the system 0 in = [0, 27r] with periodic boundary conditions and initial conditions 
parameterized by 

Mo(a;; fi) = fi | sin(2 x) | +0.1, 

where [i £ [0,1]. In this setting, the solution develops a shock that travels with the velocity a= 0.55/r. A 
dictionary T> is constructed by sampling the parameters {0, 0.2,0.4, 0.6,1.0} (r = 5) and the solution sought 
for the predictive case /i* = 0.5. A shock appears at t = 1. We display the solutions obtained by L 2 -norm, 
I/ 1 -norm by LP and IRLS and the Hubert-IRLS minimization procedures for t=j<l,t=J and t = 7r in 
Figures [TTT| and [T~T] Before the shock appears, there is almost no difference between the four solutions. The 


7T 

T = ■ 

4 




(a) Solutions 


(b) Zoom near a maximum 


T = 


7r 

4 



Figure 13: Unsteady Burgers’ equation: predicted solutions at target parameter fx* = 0.5 at t = ^ 

approach based on minimizing the L 2 -norm is even slightly better, as it can be observed from the two zooms 
in Figure [13] However, the situation after the shock appears is very different, as observed in Figure [14] the 
L 2 -norm solution is clearly the worst one with large oscillations. The L 1 -norm-type solutions are all close to 
each other and the shock is rather well reproduced with, however, an artifact that develops for longer times, 
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as seen at t = n. Nevertheless, the Zd-norm-type solutions are within the bounds of the “exact" solution, 
and no large oscillation develops. 



(a) (b) 

Figure 14: Unsteady Burgers’ equation: predicted solutions at target parameter /z* = 0.5 at t = § (left) 
and t = tt (right) 

In a second set of numerical experiments, we consider the influence of the sampling parameter set included 
in the dictionary V. We consider two dictionaries T>i = {0.4,0.45, 0.55,0.6} and V 0 = {0,0.2, 0.4,0.45,0.55, 0.6,1.0}, 
for the same target value of /z* = 0.5. These choices amounts to selecting samples close to the target value 
0.5 while varying elements of the dictionary that are not close to 0.5. We see that refining the dictionary 




Figure 15: Unsteady Burgers’ equation: predicted solutions at target parameter /z* = 0.5 at t = tt for two 
dictionaries associated with two samples of the parameter domain V 

has a positive influence as the target solution is much closer to the dictionary elements. This is confirmed 
by additional experiments where the samples of /z used to generate the dictionary were more numerous and 
closer to 0.5 (not reported here). However, keeping values of /z that are far from 0.5 has an effect, clearly 
negative for the L 2 -norm optimization. The L-norm-type solutions are however unaffected by the presence 
of these “outliers" in the dictionary, similarly as in the simple regression case reported in Section |4~T| Overall, 
the most accurate predictions are based on the minimization of the Hubert function by IRLS. 
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4.3.2 Euler equations 

The one-dimensional compressible Euler equations are considered on Q = [0,1] 



d 


pu 


+ ~- P« + P 
dx \u(E+p). 


= 0 , 


for which U = (p, pu, E) T and the pressure is given by 


P = (7 - !) ( E - 2 Pu 2 


(15a) 


(15b) 


with 7 = 1.4. 

This problem is parameterized by the initial conditions Uo(x',p). To define the parameterized initial 
conditions of the problem, the Lax and Sod cases are first introduced as follows. 

The state Us oA (x) is defined by the primal physical quantities: 


and Uh ax (x) defined by 


Lsod(^) 


1 / Lax(^') — 


p = 1 if X < 0.5, 
u = 0.0 

p = 1.0 if x < 0.5, 


p = 0.445 if x < 0.5, 
u = 0.698 if x < 0.5, 


0.125 otherwise, 
0.1 otherwise, 


0.5 otherwise, 
0.0 otherwise, 


(15c) 


(15d) 


p = 3.528 if x < 0.5, 0.571 otherwise. 


The Sod condition presents a fan, followed by a contact and a shock. For the density and the pressure, 
the solution behaves monotonically, and the contact is moderate. The Lax solution has a very different 
behavior and the contact is much stronger. This is depicted in Figure [16] where the two solutions are shown 
for t = 0.16. 

The initial condition are parameterized for p £ [0,1] as 


V 0 {x\n) = pV SoA (x) + (1 - n)V Lax {x) 


(15e) 


and the conservative initial variables Uq{x\ p) constructed from Vo(a:;p). 

In the subsequent numerical experiments, two strategies are exploited to construct, from the dictionary 
T >, the approximation u"(p) of the solution at each time step n: 

• Either we reconstruct together the discretized density vectors p, momentum m = pu and energy E, 
i.e. the state variable at time t n using only one coefficient vector a n = (a", • • • , a") 



J2 a ? u "(Pi)- 
3 =1 


(16) 


Here the {a"}j =1 are obtained by minimizing J on the density components of the state because the 
density enable to detect fans, contact discontinuities and shocks, contrarily to pressure and velocity 
which are constant across contact waves. Doing so we expect to control better the numerical oscillations, 
if any, than with the other physical variables. Similar arguments could be applied with the other 
conserved variables as well. 


• Alternatively, we reconstruct each conserved variable separately 
N N 

P" 


N 


Y, a n p P n (Pj), m" * ^ "(p 3 -), E" * £ a£E”(p,). 

3 =1 3 =1 i=1 

where the minimization procedures are done independently on each conserved variable. 


( 17 ) 
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Figure 16: One-dimensional Euler equations: density, velocity and pressure for the Lax and Sod problems 


In order to test these approaches, the PDE is discretized by finite volumes using a discretization resulting 
in Np = 3000 dofs. The parameter range V = {0.0,0.2, 0.4,0.5, 0.8,1} is considered together with a target 
H* = 0.6. The results using the first strategy, see eq. (161, are displayed in Figure 17 and those using the 
second strategy, see eq. reported in Figure [l8| 

From both figures, we can see that the overall structure of the solutions is correct. Nevertheless, there 
are differences that can be highlighted. From Figure [T7j we can observe that the density predictions, besides 
an undershoot at the shock, are good for all minimization procedures except the L 2 -norm-based one which 
is oscillatory. However, we cannot recover correct values of the initial velocity (see left boundary), because 
there is no reason to believe that the coefficient a , evaluated from the density only , will also be correct 
for the momentum. A careful observation of the pressure plot also reveals the same behavior which is not 
satisfactory. For the same reason, if any other single variable is used for a global approximation of each 
conservative variables, there no reason why better qualitative results could be obtained. 

This problem does not occur with the second strategy for the reconstruction the correct initial 

values are recovered. The four minimization procedures behave approximately the same, and we have a 
slight undershoot/overshoot at the foot of the shock. This does not appear for the L 2 -norm minimization, 
but the density is less satisfactory between the contact and the shock in that case. All methods have slight 
problems on the velocity, between the contact and the shock. Overall, the Huber function minimization 
seems to produce the most satisfactory results. 
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(a) Density p 




Figure 17: 
expansion 


One-dimensional Euler equations: 


predicted solutions with strategy (16) based on a single 


In order to obtain these results we have been faced to the following issue. Take the momentum, for 
example. For at least half of the mesh points, its value is 0, and for half of the points, its value is set to a 
constant. Hence, the matrix A used in the minimization procedure and built on the momentum dictionary 
has a rank 2 only. The same is true for the other variables, and we are looking here for r coefficients. Two 
approaches can be followed to address this issue. The first one relies on Gram-Schmidt orthogonalization 
of the solutions prior to their use as a basis for the solution. The second approach, followed here, consists 
into perturbing infinitesimally and randomly the matrices involved in the procedure, so Aij is replaced by 
Ajj +£ij. The distribution of is uniform. This has the effect of giving the maximum possible rank to the 
perturbed matrisj^] We have expressed that e l: j should depend on the variable, we have chosen 

£'lj — £ij bref 

where L re f is the difference between the minimum and the maximum, over the dictionary, of the considered 
variable. Choosing the same £i 3 is taken for all variables, this has the effect of increasing the amplitude of 
the oscillations after then shock, and our experience is that the L 2 -norm procedure is more sensitive to the 

1 We have had to use the same technique for Burgers’ equation in section 4.3.1 
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(a) Density p 




Figure 18: One-dimensional Euler equations: predicted solutions with strategy ( |17| ) based on multiple 
expansions 


choice of L re f. 

All this being said, the solution using three distinct coefficients obtained independently is of significantly 
much better quality than the one using only one expansion. 


5 Conclusion 

A novel model reduction that relies on a dictionary approach is developed and tested on several steady and 
unsteady hyperbolic problems. All of the solutions of the problem tested are parameterized and have regions 
of their spatial domain with discontinuities, leading to solutions with very distinct behaviors, such as different 
wave speeds and shock locations, making them challenging to reduce using classical projection-based model 
reduction techniques. 

To address this challenge, the proposed approach based on a dictionary of solutions is coupled with a 
functional minimization. The analysis and numerical experiments conducted in this work show that the pro¬ 
posed approach is robust (at least for one-dimensional problems) and performs best when the functional is of 
L 1 -norm-type. Among those, the Huber functional is found to be the most robust in all test cases. For all the 
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functional considered, effective numerical algorithms for their minimization are considered, both in the lin¬ 
ear and nonlinear cases. Algorithms based on convex programming are the most computationally expensive 
ones and, as a result, approaches based on iteratively solving L 2 -norm minimization problems are consid¬ 
ered. These approaches are found to be inexpensive and lead to similar solutions as convex programming. 
Furthermore, hyper-reduction approaches developed for the solution of L 2 -norm minimization 133 HIES 
can be readily applied for their efficient solution. This will be the subject of further work. Another extension 
of the present work will be the development of appropriate parameter sampling techniques together with 
error estimators to select the dictionary elements. 

For unsteady systems, in the present paper, only dictionary members computed at the same time instant 
are considered in the reduced approximation. Variants of the proposed framework will be explored. In 
particular dictionary members associated with multiple time instants will be considered as well, resulting in 
local ROM approaches El M [13- 


Acknowledgements 

The first author has been funded in part by the MECASIF project (2013-2017) funded by the French "Fonds 
Unique Interministeriel" and SNF grant //- 200021_153604 of the Swiss National Foundation. The last 
author would like to acknowledge partial support by the Army Research Laboratory through the Army High 
Performance Computing Research Center under Cooperative Agreement W911NF- 07-2-0027, and partial 
support by the Office of Naval Research under grant no. N00014-11-1-0707. This document does not 
necessarily reflect the position of these institutions, and no official endorsement should be inferred. 


References 

[1] P A LeGresley and J J Alonso. Airfoil design optimization using reduced order models based on proper 
orthogonal decomposition. AIAA Paper 2000-2545 Fluids 2000 Conference and Exhibit, Denver, CO, 
pages 1-14, 2000. 

[2] T Bui-Thanh, K Willcox, and O Ghattas. Parametric reduced-order models for probabilistic analysis 
of unsteady aerodynamic applications. AIAA Journal, 46(10):2520-2529, 2008. 

[3] D Amsallem, J Cortial, and C Farhat. Toward real-time computational-fluid-dynamics-based aeroelastic 
computations using a database of reduced-order information. AIAA Journal, 48(9):2029-2037, 2010. 

[4] D Amsallem, M J Zahr, Y Choi, and C Farhat. Design Optimization Using Hyper-Reduced-Order 
Models. Structural and Multidisciplinary Optimization, pages 1-22, 2014. 

[5] J.L. Lumley. The structure of inhomogeneous turbulent flows. In A.M. Iaglom and V.I. Tatarski, editors, 
Atmospheric turbulence and Radio wave propagation, pages 221-227, Moscow, 1967. Nauka. 

[6] L Sirovich. Turbulence and the dynamics of coherent structures. Part I: coherent structures. Quarterly 
of applied mathematics, 45(3) :561 571, 1987. 

[7] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 
1990. 

[8] Y. Maday and E. M. Rpnquist. A reduced-basis element method. J. Sci. Comput., 17(1-4):447-459, 
2002 . 

[9] K. Willcox. Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. 
Comput. Fluids, 35(2):208-226, 2006. 


25 


[10] K. Veroy and A.T. Patera. Certified real-time solution of the parametrized steady incompressible navier 
stokes equations: rigorous reduced-basis a posteriori error bounds. International Journal on Numerical 
Methods in Fluids, 47:773-788, 2005. 

[11] C.W. Rowley, T. Colonius, and R.M. Murray. Model reduction for compressible flow using POD and 
Galerkin projection. Physica D: Non linear phenomena, 189(1-2):115-129, 2004. 

[12] K. Carlberg, C. Farhat, J. Cortal, and D. Amsallem. The GNAT method for non linear model reduction: 
effective implementation and application to cpmputational fluid dynamics and turbulent flows. J. 
Comput. Phys., 242:623-647, 2013. 

[13] D Amsallem, M J Zahr, and K Washabaugh. Fast Local Reduced Basis Updates for the Efficient 
Reduction of Nonlinear Systems with Hyper-Reduction . Accepted for publication, Special issue on Model 
Reduction of Parameterized Systems (MoRePaS), Advances in Computational Mathematics, pages 1-34, 
2014. 

[14] D Amsallem and C Farhat. Interpolation method for adapting reduced-order models and application 
to aeroelasticity. AIAA Journal, 46(7):1803-1813, 2008. 

[15] M.F. Barone, D.J. Segalman, and H.K. Thornquistand I. Kalashn ikova. Galerkin reduced order mod¬ 
els for compressible flow with structural interaction. Technical report, AIAA 46th Aerospace Science 
Meeting and Exhibit, 2008. AIAA 2008-0612, Reno. 

[16] M.F. Barone, I. Kalashnikova, D.J. Segalman, and H. Thornqu ist. Stable galerkin reduced order models 
for linearized compressible gfows. J. Comput. Phys., 288:1932-1946, 2009. 

[17] M Lesoinne, M Sarkis, U Hetmaniuk, and C Farhat. A linearized method for the frequency analysis of 
three-dimensional fluid/structure interaction problems in all flow regimes. Computer Methods in Applied 
Mechanics and Engineering, 190:3121-3146, 2001. 

[18] Yvon Maday and Benjamin Stamm. Locally adaptive greedy approximations for anisotropic parameter 
reduced basis spaces. SIAM J. Sci. Comput., 35(6):2417-2441, 2013. 

[19] S Kaulmann and B Haasdonk. Online Greedy Reduced Basis Construction Using Dictionaries. VI 
International Conference on Adaptive Modelling and Simulation (ADMOS 2013), 2013. Moitinho de 
Almeida, Jose Paulo Baptista and Diez, Pedro and Tiago, Carlos and Pares, Nuria (Eds.),. 

[20] C. Prud’homme, D.V. Rovas, K. Veroy, L. Machiels, Y. Maday, A.T. Patera, and G. Turicini. Reliable 
real-time solution of parametrised partial differential equations: reduced basis output bound methods. 
Journal of Fluids Engineering, 124:70-80, 2002. 

[21] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera. An "empirical interpolation" method: 
Application to efficient reduced-basis discretization of partial differential equations. C. R., Math., 
Acad. Sci. Paris, 339(9):667-672, 2004. 

[22] A. Buffa, Y. Maday, A. T. Patera, C. Prud’homme, and G. Turinici. A priori convergence of the 
greedy algorithm for the parametrized reduced basis method. ESAIM, Math. Model. Numer. Anal., 
46(3):595-603, 2012. 

[23] Y. Chen, J. S. Hesthaven, Y. Maday, and J. Rodriguez. A monotonic evaluation of lower bounds for 
inf-sup stability constants in the frame of reduced basis approximations. C. R., Math., Acad. Sci. Paris, 
346(23-24): 1295-1300, 2008. 

[24] A Paul-Dubois-Taine and D Amsallem. An adaptive and efficient greedy procedure for the optimal train¬ 
ing of parametric reduced-order models. International Journal for Numerical Methods in Engineering, 
pages 1-31, September 2014. 


26 



[25] D. Ryckelynck. A priori hyperreduction: an adaptive approach. J. Comput. Phys., 202:346-366, 2005. 

[26] J.L. Guermond and B. Popov. lA-approximation of stationnary Hamilton-Jacobi equations. SIAM J. 
Numer. Anal., 47(l):339-362, 2008. 

[27] J.L. Guermond, F. Marpeau, and B. Popov. A fast algorithm for solving first-order PDEs by L 1 
minimization. Communications in mathematical Sciences, 6(1):199-216, 2008. 

[28] C. T. Kelley and D. E. Keyes. Convergence Analysis of Pseudo-Transient Continuation. SIAM Journal 
of Numerical Analysis, 35(2):508-523, 1998. 

[29] Ingrid Daubechies, Ronald DeVore, Massimo Fornasier, and C Sinan Gunturk. Iteratively re-weighted 
least squares minimization for sparse recovery, http://arxiv.org/abs/0807.0575 July 2008. 

[30] Peter J Huber and Elvezio M Ronchetti. Robust Statistics. John Wiley & Sons, September 2011. 

[31] S Chaturantabut and DC Sorensen. Nonlinear model reduction via discrete empirical interpolation. 
SIAM Journal on Scientific Computing, 32(5):2737-2764, 2010. 

[32] K Carlberg, C Bou-Mosleh, and C Farhat. Efficient non-linear model reduction via a least-squares 
Petrov-Galerkin projection and compressive tensor approximations. International Journal for Numerical 
Methods in Engineering, 86(2):155-181, 2011. 

[33] D Amsallem, M J Zahr, and C Farhat. Nonlinear model order reduction based on local reduced-order 
bases. International Journal for Numerical Methods in Engineering, 92(10):891 916, 2012. 

[34] M Dihlmann, M Drohmann, and B Haasdonk. Model reduction of parametrized evolution problems 
using the reduced basis method with adaptive time-partitioning. Proc. of ADMOS, 2011, 2011. 

[35] Jorge Nocedal and S J Wright. Numerical optimization. Springer, December 2006. 

A Algorithms 

This section reviews the minimization algorithms that have been used in this study. 

A.l L 2 -norm minimization 

A. 1.1 Linear Case 

The least-squares problem is defined for a skinny matrix A = MD £ M. Nxr and a vector b £ R w as 

min || Az + b|||. (18) 

Z 

One possible solution is by QR decomposition as described in Algorithm [l] 


Algorithm 1 Linear L 2 -norm minimization (Least-squares) by QR decomposition 

Input: Matrix A and vector b 
Output: Solution z 
1: Compute the QR decomposition of A 

A = QR 


2: Let z = Q 1 R^b 
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A.1.2 Nonlinear Case 


This algorithm is used in the case drawn in Remark |3.2| The following nonlinear least-squares problem can 
be solved by the Gauss-Newton method as described in Algorithm [2] 

min ||r(z)|||. (19) 

Z 

Another approach to solve this problem is by the Levenberg-Marquardt procedure [35 . 


Algorithm 2 Nonlinear L 2 -norm minimization by the Gauss-Newton method 

Input: Residual function r(-) and associated Jacobian J(-), dictionary D. initial guess z°, tolerance for 
convergence e 
Output: Solution z 


1 : l = 0 

2: Compute r° = r(Dz°) and W° = J(Dz°)D 
3 : while (W l ) T r l >e 


(W°) T r° 


do 


4 : Solve the linear L 2 -norm minimization problem using Algorithm [l] with arguments ~W l and 


Az l = argmin || W*y + r l 
y 


2 

2 


5 : = z} + Az ( 

6: Compute r /+1 = r(Dz i+1 ) and W ,+1 = J(Dz i+1 )D 

7 : 1 = 1 + 1 

8: end while 

9 : Z = Z l 


A.2 L^-norm minimization 

A.2.1 Linear Case 

The linear L 1 -norm minimization problem is defined for a skinny matrix A = MD £ R Arxr and a vector 
b e as 

min |j Az + b|| ]_. (20) 

Z 

This problem can be recast as a Linear Program as described in Algorithm [3j 


Algorithm 3 Linear L 1 -norm minimization by Linear Programming 

Input: Matrix A and vector b 
Output: Solution z 
1: Solve the linear program 

(z*, s*,t*) = argmin l T (s + t) 

Z,S,t 

s.t. Az — s + t = b 
s > 0 

t > 0 


2: Let z = z* 


An issue associated with the solution of (20) is the fact that there are 2 N+r variables and 3 N constraints, 
among which N are equality constraints. N is the number of degrees of freedom in the high-dimensional 
problem and can be very large for fine discretization problems, rendering the linear program solution in¬ 
tractable. 
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Alternatively, the L 1 -norm minimization problem (201 can be solved by Iteratively Reweighted Least- 
Squares (IRLS) [ 22 ] • This approach proceeds iteratively by solving a sequence of weighted least-squares 
problem. An advantage of this approach is that its implementation can rely entirely on existing least- 
squares solvers. Furthermore, its complexity is similar to that of the L 2 -norm minimization problem. The 
procedure is presented in Algorithm [4j 


Algorithm 4 Linear L 1 -norm minimization by Iteratively Reweighted Least-Squares (IRLS) 

Input: Matrix A and vector b, initial guess z° 

Output: Solution z 
1 : 1 = 0 

2: Compute r° = ADz° + b and W° = AD 
3: while l = 0 OR ||Az i_1 ||i > e(l + ||z i_1 ||i) do 

4: Compute the weights 7} = diag ^ | r\ | 

5: Solve the linear L 2 minimization problem using Algorithm [l] with arguments 7 l ~W l and 7 l r l 

Az l =argmin||Z i W'y + ZV || 2 
y 


6 : Z l+1 = Z l + Az l 

7: Compute r i+1 = ADz i+1 + b and W 1+1 = AD 

8 : 1 = 1 + 1 

9: end while 

10: Z = Z l 


A.2.2 Nonlinear Case 

The following nonlinear iA-norm minimization problem can be solved by a Gauss-Newton-like procedure. 

min ll r ( z )lli- ( 21 ) 

Z 

The approach relies on the solution of a sequence of linear LCnorm minimization problems. Algorithm [ 5 ] 
describes the approach when it relies on Linear Programming and Algorithm [6] when it relies on IRLS. 


Algorithm 5 Nonlinear Zd-norm minimization by the Gauss-Newton method with LP 

Input: Residual function r(-) and associated Jacobian J(-), dictionary D, initial guess z°, tolerance for 
convergence e 
Output: Solution z 
1 : l = 0 

2: Compute r° = r(Dz°) and W° = J(Dz°)D 
3: while 1 = 0 or |||W'Az' + r z ||i - ||r z ||i| > eHrl! do 

4: Solve the linear ZJ-norm minimization problem using Algorithm [ 3 ] with arguments ~W l and R 

Az l = argmin ||W*y + r*||i 
y 


5. Z / + l 2} _|_ /\y} 

6: Compute r /+1 = r(Dz i+1 ) and W l+1 = J(Dz z+1 )D 

7: 1 = 1 + 1 

8: end while 

9: Z = Z l 
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Algorithm 6 Nonlinear L [ -norm minimization by the Gauss-Newton method with Iteratively Reweighted 
Least-Squares 

Input: Residual function r(-) and associated Jacobian J(-), dictionary D, initial guess z°, tolerance for 
convergence e 
Output: Solution z 
1 : l = 0 

2: Compute r° = r(Dz°) and W° = J(Dz°)D 
3: while l = 0 OR |jAz i_1 |i > e(l + || z z “ 1 1|l) do 
4: Compute the weights Z l = diag ^|r-| 

5: Solve the linear L 2 -norm minimization problem using Algorithm [l] with arguments Z l 'W l and Z l r l 

Az l = argmin||Z i W i y + ZV||^ 
y 

6 : Z l + 1 =Z l + Az l 

7: Compute v l+1 = r(Dz i+1 ) and W ,+1 = J(Dz /+1 )D 

8 : 1 = 1 + 1 

9: end while 

10: Z = Z l 


A.3 Huber function minimization 

An issue with L 1 -norm minimization is the fact that the function z i —> ||z||i is non differentiable at z = 0, 
causing potential difficulties in the numerical solution of the minimizer, as shown in the numerical results 
of Section [4. 2. 2| Alternatively, the minimization of the Huber function can be used. This function behaves 
similarly to x 2 for small values of x and as |x| for large values of x. It is also differentiable everywhere. The 
Huber function ipM is defined as: 


4>m{x) 


x 2 if |x| < M 

M(|x| — M ) otherwise, 


The L '-norm minimization problem is then replaced by 


( 22 ) 


N 

min V'^ M (r i (z)) (23) 

Z Z ' 

2—1 

This problem can also be solved by the IRLS approach as described in Algorithm [ 7 ] This approach requires 
choosing an appropriate value for M. In the present work, the following choice has been found to be robust 
across all applications 

M = e 2 max(l, max(|ri|)) 

with e -2 = 10 -6 . 

The L 2 -norm, L [ -norm and Huber function minimizations can be all recast with different choices of 
functions (j> as 

N 

( 24 ) 

z z —' 

2=1 

Figure [T9] compares the different functions (j) involved. 
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Algorithm 7 Nonlinear Huber function minimization by Iteratively Reweighted Least-Squares 

Input: Residual function r(-) and associated Jacobian J(-), dictionary D. initial guess z°, tolerance for 
convergence e 
Output: Solution z 
1 : l = 0 

2: Compute r° = r(Dz°) and W° = J(Dz°)D 
3: while l = 0 OR ||Az i_1 |i > e(l + || z z ~ 1 1|l) do 

4: Compute the weights Z l = diag < M) + M\r\\ 2 S(ri > M)j 

5: Let M = €2 max(l,max(|rj|)) 

6: Solve the linear L 2 -norm minimization problem using Algorithm [l] and arguments Z l ~W l and Z l r l 

Aq' =argmin||Z i W'y + ZV ||2 
y 


7 : q' +1 = q z + Aq' 

8: Compute v l+1 = r(Dz I+1 ) and W ,+1 = J(Dz /+1 )D 

9: 1 = 1 + 1 

10: end while 

11: Z = Z l 



X 


Figure 19: Comparison of the L 2 , L 1 and Huber function norms 
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